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Abstract —We present a new method for time delay estimation 
using band limited frequency domain data representing the port 
responses of interconnect structures. The approach is based on 
the recently developed by the authors spectrally accurate method 
for causality characterization that employs SVD-based causal 
Fourier continuations. The time delay extraction is constructed 
by incorporating a linearly varying phase factor to the system 
of equations that determines Fourier coefficients. The method 
is capable of determining time delay using data affected by 
noise or approximation errors that come from measurements 
or numerical simulations. It can also be employed when only a 
limited number of frequency responses is available. The technique 
can be extended to multi-port and mixed mode networks. Several 
analytical and simulated examples are used to demonstrate the 
accuracy and strength of the proposed technique. 

Index Terms —time delay, delay estimation, causality, disper¬ 
sion relations, singular value decomposition, SVD-based causal 
Fourier continuation, high speed interconnects. 


I. Introduction 

Identification and extraction of time delay is an important 
research problem in signal processing and has applications in 
many fields including radar 1251 . sonar |33) . |f8l , ultrasonics 
Q0), microwave imaging (26), geophysics l22l . seismology 
(39), 129) . wireless communications |36| as well as model¬ 
ing of passive structures in electronic systems, in particular, 
transmission line modeling m no, transient simulation 
of interconnects ll24l and co-simulation of passive structures 
with active devices in a time domain using SPICE. Passive 
structures in electronic systems have been traditionally an¬ 
alyzed in the frequency domain, while transient simulations 
are performed in the time domain using suitable models that 
accurately capture the relevant electromagnetic phenomena. 
The models are obtained from either direct measurements or 
electromagnetic simulations. Interconnect models are typically 
approximated by rational transfer functions using the vector 
fitting algorithm in various implementations ED, Q2), Go), 
E2. CD, E), 0, which is the standard macromodeling 
approach. As clock frequencies increase, the size of passive 
structures becomes of the same order as the signal wavelength 
at the operating frequency, which causes the distributed effects 
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such as time delay to play a significant role in the time domain 
simulations. For this reason, time delay has to be included in 
macromodeling, in particular, when causality is analyzed. The 
connection between causality and time delay is in the fact 
that time delays can pull a non-causal signal into the causal 
region or vice versa pull a causal signal into the non-causal 
region, while causality, in turn, can be expressed in terms of 
the Hilbert transform ED, m, EQ. Several approaches can 
be used to extract delays in the frequency domain, for example, 
using the Hilbert transform 01, G2), GD, the minimum 
phase all-pass decomposition (28), GZ), (24), incorporating an 
optimal time delay into the vector fitting algorithm [18), [11], 
employing a modified Lie approximation to develop a passive 
and compact macromodel GO), using a Gabor transform 
to develop delayed rational function macromodels for long 
interconnects da, 0 or conducting a probabilistic analysis 
of the cepstrum in the presence of noise ED. In the time 
domain, delayed rational functions 0, G) can be employed 
to extract delays. In this paper, a novel approach is proposed in 
which time delay is determined in the frequency domain using 
a causality argument. Causality is verified using the SVD- 
based causal Fourier continuation method developed by the 
authors 0. 0, while the time delay presence is incorporated 
by a linearly varying phase factor to the system of equations 
that determines Fourier coefficients. Preliminary results are 
reported in |5) . 


The rest of the paper is organized as follows. Section 071 
provides a background on causality for linear time-translation 
invariant systems and dispersion relations. In Section [HI] 
we show main steps in the derivation of causal Fourier 
continuations using truncated singular value decomposition 
(SVD) method that was developed to access causality. We 
also provide error estimates that take into account a possible 
presence of noise in data. Section m extends the causal¬ 
ity characterization method to develop a technique for time 
delay extraction. The proposed method is tested in Section 
m using several analytic and simulated examples. We also 
analyze the performance of the algorithm when only a lim¬ 
ited number of frequency responses is available and when 
noise/approximation errors are present in data. In Section [VT] 
we present our conclusions. The Appendix section is devoted 
to formulation of error bounds for the causality characteriza¬ 
tion method based on causal Fourier continuations. 
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II. Causality of Linear Time-Invariant Systems 


Consider a linear and time-invariant physical system with 
the impulse response h(t) subject to a time-dependent input 
fit), to which it responds by an output x(t). Denote by 

/ OO 

h{r) e~ iwT dr (1) 

-OO 


the Fourier transform of h{t), which is also called the transfer 
function. 

The system is causal if the output cannot precede the input, 
i.e. if f{t) = 0 for t < T, the same must be true for xit). 
This primitive causality condition in the time domain implies 
hit) = 0, t < 0. Hence, domain of integration in ijT} can be 
reduced to [0, oo). 

Assume H{w) £ L 2 (K). Then starting from Cauchy’s 
theorem and using contour integration, one can show [311 that 
for any point w on the real axis, H[w) can be writterjj as 


where 



( 2 ) 


denotes Cauchy’s principal value. Separating the real and 
imaginary parts of (0, we get 


Re 



Im H(w) . 

- aw 

w — w' 


(3) 


Im H (w) 



ReHiw) , 
- —aw 

w — w' 


(4) 


Expressions ([3} and © are called the dispersion relations or 
Kramers-Kronig relations. They show that R eH and Im // 
are not independent functions, but instead they are related 
to each other: R eH at one frequency depends on Im H at 
all frequencies, and vice versa. This implies that if one of 
the functions Re H or Im H is square integrable and known, 
then the other one can be completely determined by causality. 
Recalling the definition of the Hilbert transform. 


i r°° 

H[u(w)] = -/ 

^ J — oo 


u(w') 
W — w' 


dw ', 


we see that R eH and I in H are Hilbert transforms of each 
other, i.e. 


ReHiw) = r H[RnH (w)], ImiT(m) = —H [Re H(w)\ ■ 

(5) 

In other words. Re H or Im H form a Hilbert transform pair. 
Dispersion relations provide the causality condition in the 
frequency domain. 

Evaluation of the Hilbert transform requires integration 
on (—oo,oo), which can be reduced to [0, oo) by spectrum 
symmetry of Hiw) if h(t) is real valued. In practice, only 
a limited number of discrete values of Hiw) is available 
on [ WminjWmax ]■ Thus, the domain of integration has to be 
truncated. This usually causes serious boundary artifacts due to 


'Please note that we use an opposite sign of the exponent in the definition 
of the Fourier transform than in ED. 


the lack of out-of-band frequency responses. To reduce or even 
completely remove boundary artifacts, the authors recently 
developed periodic polynomial ID, 0 and causal Fourier 
continuation ii, a, respectively, based methods for causality 
characterization. The approach was motivated by the example 
Hiw) = e~ law , a > 0, that is not square integrable but still 
satisfies the dispersion relations. The causality characterization 
method based on causal Fourier continuations allows one to 
construct highly accurate approximations of a given transfer 
function on the original frequency interval [w mm ,t« m oi] with 
the uniform error that decreases as the number of Fourier 
coefficients increases. The technique is applicable to both 
baseband and bandpass cases and capable of detecting very 
small localized causality violations. The method can also be 
extended to multidimensional cases. 

In the next section, for completeness of presentation, we 
show main steps in the derivation of the causal Fourier 
continuation method that can be used to access causality of 
a given transfer function whose values are available at a 
discrete set of frequencies. We also provide upper bounds 
of reconstruction error between the given function and its 
causal Fourier continuation. We use these error estimates to 
understand how to extract time delay when data with different 
resolutions are available and when data are affected by noise 
or other approximation errors. 

III. Causal Fourier Continuations 

Consider a transfer function H (u>) = ReH+i Im H , whose 
N discrete values are available on [w rn i n , w max ], w m i n > 0. 
For real-valued impulse response functions hit), ReH and 
ImH are even and odd functions, respectively. This implies 
that Hiw) has values on [— w max , — w m in\ by spectrum sym¬ 
metry. For convenience, we rescale the frequency interval 
[-W max ,w max \ to [-0.5, 0.5] by the substitution x = x^-w, 
so the rescaled transfer function H(x) is defined on the unit 
length interval with N values where N = 2 N — 1 or TV = 2 N 
depending if H(x) is available at x = 0 or not. Both baseband 
and bandpass cases can be considered. 

The idea of a causal Fourier continuation is to construct an 
accurate Fourier series approximation of IT(x) by allowing 
the Fourier series to be periodic and causal in an extended 
domain. The result is the Fourier continuation of H that we 
denote by C(7T), and it is defined by 

M 

CiH)ix)= J2 ^e~^ kx , (6) 

k=—M+l 

for even number 2 M of terms, whereas for odd number 2M+1 
of terms, the index k varies from —M to M. Throughout this 
paper, we assume that the number M of Fourier coefficients is 
even, for simplicity. When M is odd, analogous results can be 
formulated. Here b is the period of approximation. For SVD- 
based periodic continuations b is normally chosen as twice the 
length of the domain on which function H is given though the 
value b = 2 is not necessarily optimal. The optimal value b 
depends on a function being approximated. In practice, several 
values b £ (1,4) may be tried to get a better reconstruction 
of H{x) with a Fourier series. 
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Functions (}>k(x) = e~^ kx , fc e Z, form a complete or¬ 
thogonal basis in L 2 [— |, |]. It can be shown that 'H{<j> fc (x)} = 
isgn(k)<j)k(x), which implies that functions {4>k( x)} are the 
eigenfunctions of the Hilbert transform H with associated 
eigenvalues ±i with x € [— |, |]. For a causal periodic 
continuation, according to (0, we need Im C(H)(x) to be the 
Hilbert transform of — ReC(fT)(x). It can be shown 0 that 
this implies ak = 0 for k < 0 in 0. Hence, a causal Fourier 
continuation has the form 

M 

C(H)(x) = J2 a k (t>k{x). (7) 

fe=l 

Evaluating H(x) at points x ? , j = 1,..., TV, Xj £ [—0.5,0.5], 
produces a complex valued system 

M 

T: oik<t>k(xj) = H(xj) (8) 

fc=i 

with TV equations for M unknowns o^, fc = 1,..., M, TV > 
M. If TV > M, the system ® is overdetermined and has to 
be solved in the least squares sense. When Fourier coefficients 
at are computed, formula © provides reconstruction of H(x) 
on [—0.5,0.5]. The least squares problem is extremely ill- 
conditioned. However, it can be regularized using a truncated 
SVD method when singular values below some cutoff toler¬ 
ance £ close to the machine precision are being discarded. To 
have a better control on ill-conditioning of matrix problem ©, 
more data points TV than the Fourier coefficients M should 
be used. We use at least TV = 2M as an effective way to 
obtain an accurate and reliable approximation of H(x) over 
the interval [—0.5, 0.5]. This relation corresponds!! to TV = M, 
where TV is the number of data points available originally on 

IXmini 'OJmax]- 

Since Re H (x) and Im H (x) are even and odd functions of 
x, respectively, the Fourier coefficients 

1 f b / 2 - 

a k = -r / H(x)(j>k(x)dx, k = 1,... ,M, 

b J-b/2 

are real. Here " denotes the complex conjugate. To ensure that 
numerically computed Fourier coefficients ak are real, instead 
of solving complex-valued system @, one can separate the 
real and imaginary parts of C(H)(xj ) to obtain real-valued 
system 

M 

R eC(ff)(*i) = £ ak Re^fe(xj), 

fe=t (9) 

M v ’ 

Im C(H)(xj) = Im (j>k{xj). 

k =1 

We show in (4) that real formulation © provides slightly more 
accurate results than complex. 

To access the quality of approximation of H(x ) with its 
causal Fourier continuation C(II) (x), we introduce reconstruc¬ 
tion errors Er(x) and Ej(x), 

E r {x) = ReiT(x) - ReC(TT)(x), (10) 

“In 31, N denoted the number of points on [—0.5, 0.5], while in this work 
N is the number of points on [0, 0.5] or originally on [w m i n , 


Ej(x) = ImfT(x) — ImC(iT)(x) (11) 

on the original interval [—0.5,0.5]. 

The error analysis performed in 0 (see also Appendix) 
shows that the error between H(x) and its causal Fourier 
continuation C (H + e) under the presence of a noise e, has 
the following upper bound: 

\\H — C(H + e)||z, 2 (f2) < e_F + £n + er- (12) 

Here 

e F = (1 + A 2V / 2TV(M - K))\\H - H M \\ Lao (n) (13) 

is the error due to approximation of H with a causal Fourier 
series and it decays as 0(M~ k+1 ), where k is the smoothness 
order of the transfer function II (x). 

e r = Aiy / R7^l|T?M||L 00 (n“) (14) 

is the error due to the truncation of singular values and 
it is typically small and close to the cut-off value £. As 
d indicates, ep depends on b and the function iT being 
approximated. 

e n = (1 + A 2y /2N(M - /0)||e|| Loo(f2) (15) 

is the error due to the presence of a noise or approximation 
errors in the given data and it shows a level of causality 
violation. In practice the size of e n is close to the size of noise 
in data. Function Hm and constants Ay, A 2 and I\ are defined 
in Appendix. These constants depend only on the continuation 
parameters TV, M, b and £ as well as location of discrete points 
Xj, and not on the function H. 

The error bound (IT2l) shows that the reconstruction errors 
Er and Ei decrease as M increases due to the causal Fourier 
series approximation error with the error bound ep until either 
the level e of a noise or level ep due to truncation of singular 
values is reached. If only round-off errors are present in data, 
the errors will level off at ep. If reconstruction errors level 
off at some value e > ep as the resolution increases, the 
data are declared non-causal with the error approximately at 
the order of e. More information about the error analysis for 
the causality characterization methods based on causal Fourier 
continuations can be found in 0. 

IV. Time Delay Estimation 

The above approach for causality assessment can be trans¬ 
formed into a delay estimation algorithm by observing the 
following. Suppose that h(t) is non-zero only from time 
To > 0, and we would like to identify the time delay l\ s . 
Consider the Fourier transform H(w) of h(t): 

pOO /* OO 

H{w) = / e~ iwt h(t)dt = / e _i “* h(t)dt 

J t=To J t=To 

where we used the substitution x = aw, a = 0,5 . Introduc- 

'Wmax 

ing r = we can write 

poo 

H{w) = a / e~ lXT h{ar)dr, 

Jlii 
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or with u = t — —, we obtain 

a 7 

t r°° t 

H(w ) =ae~ ix ^ / e" ia:t /i(T 0 + au)du = ae-“^ G(®), 
Jo 

where G(ir) is the Fourier transform of a causal function with 
no time delay. This implies that when 0 < T < Tq, the transfer 
function H{x)e lx ~ is causal, but when T > To, the transfer 
function H(x) e lx ~ has a non-causal component. Therefore, 
To = To/a is the time delay for H(x), and the delay To for 
the original function H(w) is recovered by multiplying To by 
a. Since one can add any integral multiple of 2tt to xT/a, it 
is enough to restrict our investigations to the interval 


T 

0 < — < T, 
a 


27r 


Kmax 



Then for each potential time delay 0 < — < T max , we solve 
the following modified system 

M 

^Zotk&kixj) = H(xj), j = 1,... ,N (16) 

fe=l 


or its equivalent real-valued formulation. For T < To, the 
reconstruction errors Er, hj should be small and approxi¬ 
mately of the same order. As T increases and becomes greater 
than some critical transition time close to the time delay To, 
the reconstruction errors should start increase. The goal is 
to approximate T 0 . The difficulty is that the reconstruction 
errors grow gradually as T > T 0 , so transition is not sharp. 
Moreover, the order of reconstruction errors for T < T 0 
depends on the resolution of data and threshold £ used in 
the truncated SVD method, which, in turn, affects a transition 
time. In addition, a noise in data, if present, also affects when 
reconstruction errors start growing. A similar approach was 
used in l23l to estimate the time delay for square integrable 
transfer functions. In this contribution, we extend the approach 
to more general transfer functions. In addition, we use a 
different causality measure than in ll23l and take into account 
different resolutions of given data and a possible presence of 
noise. The approach can be extended to multi-port and mixed 
mode networks by applying it to each element of the transfer 
matrix. 


V. Numerical Examples 

In this section, we apply a proposed technique to several 
analytic and simulated examples when the time delay is either 
known exactly or can be estimated using other techniques. We 
also consider the effect of noise presence on the accuracy of 
timed delay estimation. 


A. Four-Pole Example 

Consider a transfer function with four poles and time delay 
To, defined by 

H(w) = e~ iwT ° H{w) (17) 


with 


H(w) = 


r 1 


r\ 


ri 


ri 


IW +P 1 IW + Pi IW + P2 IW + P2 


where n = l + 2i,p\ = l + 3i, V 2 = § + \i, P 2 = \ + 5i, and 
To = 0.25. Since the poles of H(w) are located in the upper 
half w-plane at ±3 + i and ±5 + i i , this function is causal 
as a sum of four causal transforms, and has no time delay. 
Therefore, the function H(w) is a causal function delayed with 
offset Tq. H is sampled on [0, w ma x} at N frequency points 
varying from 50 to 1500 with w rnax = 6. The real and imag- 



Fig. 1. R eH(w) and lmH(w) in the four-pole example. 


inary parts of H(w ) are shown in Fig. [I] After rescaling with 
x= T ^w and reflecting to negative frequencies, we obtain 
a rescaled transfer function H(x) defined on x G [—0.5, 0.5], 
for which we construct a causal Fourier continuation C(H) 
defined in (0 using M = N Fourier coefficients. Hence, 
the number M of Fourier coefficients also varies between 
50 and 1500. ReiT(x) and ImiT(a;) of the rescaled and 
reflected H (x) together with their causal Fourier continuations 
with M = 300 are depicted in Fig. [2] Even though given 
H (x) and its causal Fourier continuation approximation look 
indistinguishable, the actual reconstruction errors Er and Ej 
in both real and imaginary parts, that are defined in (ITOb . ( ITTb . 
are on the order of 10 -6 and they decreases as M increases 
(with M = N). For example, with M = 800, the errors are 
on the order of 10 -13 . Since both errors Er and Ej are of 
the same order, it is enough to analyze one of the errors, for 
example, Er. The results using Ei are similar. 

To estimate the time delay, we analyze the evolution of 
the ||Tij||oo, shown in Fig. 0 for various values M. Since 
the error due to a causal Fourier series approximation de¬ 
creases with M (see error bound (Il3l) ). the reconstruction 
error between the given transfer function El and its causal 
Fourier continuation C(H) also decreases as M increases until 
it either reaches the level £ of filtering of singular values or 
a level e of noise/causality violations (see error bounds (IT4l) 
and G3. respectively). For each fixed M, as time T increases, 
the errors Er and Ej first are small and about of the same 
order until some transition time close to the time delay T 0 
is approached. After that the errors grow approximately as 
a power function on the loglog scale. For smaller T, the 
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Fig. 2. R eH(x) and Im H(x) of the rescaled and reflected by symmetry 
transfer function H(x) in the four-pole example together with their causal 
Fourier continuations R eC(H) and Im C(H), respectively, with M = 300 
Fourier coefficients. 



Fig. 3. Evolution of the reconstruction error ||f?i{||oo as a function of T 
in the four-pole example. The dashed line corresponds to the exact delay 
T 0 = 0.25. 


errors are dominated by a causal Fourier series approximation 
error and then for T greater than some transition time - by 
causality violations since this value T provides large enough 
negative time delay and shifts a causal function into a non- 
causal area. A transition value T = T c , we call it a critical 
time, from a plateau region to a growth region, is different 
for each M and it decreases as the resolution or number of 
Fourier coefficients increases if the error is dominated by the 
causal Fourier series approximation error. The critical times 
T c approach the time delay To as M increases. The goal 
is to estimate To using the error curves shown in Fig. [3] 
Analyzing graphs of the error curves for M > 800, we 
observe some non-monotonic behavior at T close to To. This 
behavior is due to the filtering of the singular values below 
the threshold £ = 10” 13 that we used in our experiments. By 
increasing the value of £, the non-monotonic behavior will be 
present at smaller values of M. This suggests that portions 
of error curves close to threshold £ are affected by filtering 
and may be inaccurate and difficult to use for time delay 
estimation as we find in our experiments. To estimate critical 
times T c of transition from the plateau region to the growth 
region, we approximate the growing region by a quadratic 


function on the loglog scale. Specifically, we assume that 
InT « a 2 (ln||T_ R || 00 ) 2 +ai lnH^Hoo+oo = /(In \\E R \ loo), 
where coefficients a o, ai, and a 2 are determined in the least 
squares sense. The resulting quadratic function /(In H-TrIIoo) 
is then evaluated at the value of ||.Er||oo at T = 0 that is 
assumed to be the “most causal” time. By taking exponential 
function of the result, we find a critical transition time T c 
for a given M. This procedure produces estimates of the 
time delay Tq for various values of M. The graph of the 
critical transition times T c as a function of M is shown in 
Fig. a One can clearly see that the critical times approach 
the exact time delay To = 0.25 as M increases. A good 
approximation of Tq is achieved at M = 800. The values 



Fig. 4. Critical transition times T c in the four-pole example that approach To 
as M increases. The dashed line corresponds to the exact delay Tq = 0.25. 


of T c for M > 200 are presented in Table Q] The results 
indicate that the approximations become more accurate as M 
increases. The error with M = 900 is less than 1 %. At the 
same time, the error with M = 1500 is about 3%, which is 
due to the fact that the results in this case are more affected 
by the filtering of singular values. In the cases when M is 
high and the resulting error is not flat for T < To, instead of 
evaluating a fitted quadratic curve at the value of 11 -E/i 11 oo at 
T = 0 we evaluate it at £, the threshold of filtering singular 
values, to avoid using results affected by filtering. 


M 

T c 

M 

T c 

200 

1.4604 

700 

0.3394 

300 

1.1294 

800 

0.2529 

400 

0.9077 

900 

0.2497 

500 

0.6655 

1000 

0.2472 

600 

0.4759 

1500 

0.2576 


TABLE I 

Critical transition times T c in the four-pole example that 
APPROACH Tq = 0.25 AS M INCREASES. 


In practice the number N of samples of the transfer function 
H(w) is usually limited, which sets the bound for number 
M = N of Fourier coefficients, so it may not always be 
possible to use large enough M to obtain critical time T c 
close enough to the actual time delay Tq. A good method 
should be capable of producing an accurate approximation 
of To even with a small number of data points. We achieve 
this by employing another approach for time delay estimation. 
Using the obtained fitted quadratic error curves, we extrapolate 
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them to the value £ of filtering of singular values, which is 
typically chosen to be close to the machine precision. This 
corresponds to finding time T at which the error reaches the 
value £. This choice is natural since the errors below £ are most 
likely affected by filtering and may not be accurate enough 
to use. The results of such extrapolation are shown in Fig. 
[5] for M = 200, 400, 600 and 800. An intersection of the 
extrapolated curve corresponding to M = 200 is at a value 
T = 0.45451, which is a bit far from the exact T 0 = 0.25. At 
the same time, intersections of extrapolated curves with higher 
values of M are much closer - see Table Ullfor details. Results 



Fig. 5. HEflllcc in the four-pole example with M = 200, 400, 600 and 800 
together with their extrapolated quadratic fits. Vertical dashed line indicates 
the exact time delay To = 0.25, while horizontal dashed line indicates the 
level of filtering of singular values given by £ = 10 -13 . 


M 

To approximation 

M 

To approximation 

200 

0.45451 

700 

0.24734 

300 

0.27235 

800 

0.24969 

400 

0.25297 

900 

0.24974 

500 

0.25053 

1000 

0.24724 

600 

0.24633 

1500 

0.25759 


TABLE II 

Approximations of T 0 = 0.25 in the four-pole example using 

EXTRAPOLATION. 


shown in this table indicate that as M increases, extrapolated 
quadratic curve intersect the horizontal line the value £ at times 
closer to To. Obtained approximations of To can be averaged 
producing To « 0.24805. The approach with extrapolation 
provides a faster convergence and good approximations of To 
even for small values of M, i.e. less data points are needed to 
approximate To. 

We also consider the effect of noise on the time delay 
estimation. To study this, we impose a sine perturbation 

asin(107nr) (18) 

of various amplitudes a that we add to RefT, while keeping 
ImH unchanged. We choose N = 800 and vary a from 10 -10 
to 10 3 . The reconstruction error Tr with no perturbation for 
early times T < To is of the order of 10 -13 , as shown in 
Fig. [6] that corresponds to the level of filtering of singular 
values. When the perturbation is added, the reconstruction 
errors for T < To are higher and approximately of the order 
of a. Once some critical transition time greater than Tq is 


passed, reconstruction errors start growing and they grow at 
the same rate and coincide almost perfectly with each other. 
This observation suggests that the proposed approach can also 
be used in the cases when data have a noise, which is typical 
in real-life applications, when data have either measurement 
or simulation errors. For noise with a smaller amplitude, the 



Fig. 6. Evolution of ||Er||oo in the four-pole example with added sine 
perturbation asin(107rai). The dashed line corresponds to the exact delay 
T 0 = 0.25. 

region close to To will be less affected by noise and a bigger 
growing region will be available for fitting, so we expect better 
accuracy of time delay estimation in such cases. When more 
noise in data is present, less growing region will be available 
for fitting and extrapolation of fitted quadratic error curves 
may be less accurate. We demonstrate this by considering two 
cases: with a = 10~ 5 (noisier case) and a = 10~ 8 (less noisy 
case). The error curves with a higher amplitude a = 10 -5 are 



Fig. 7. Evolution of ||Tr||oo in the four-pole example with the added 
perturbation 10 —5 sin(107nc). The dashed line corresponds to the exact delay 
T 0 = 0.25. 

presented in Fig. [7] It is clear that the error does not become 
smaller than 10 -5 as M > 300 gets larger because of the 
noise. We use available growing regions and extrapolate fitted 
error curves to find their intersection with the horizontal line 
with value £. This gives us time T when the error reaches 
the value £ for each considered M. The results of such 
extrapolation for M = 200, 400, 600 and 800 are shown in 
Fig-® Clearly, extrapolated error curves reach value £ at times 
around T 0 but not close enough to T 0 and without established 
convergence but rather in a spread-out manner around Tq. 
Approximations of To for values of M that we investigated 
are shown in Table [III] Averaging these approximations we 
obtain = 0.26586. The extrapolated curves can be 
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Fig. 8. Error curves ||£'^|| 00 in the four-pole example with the added 
perturbation 10 -5 sin(107rx) and M = 200, 400, 600 and 800 together with 
their extrapolated quadratic fits. Vertical dashed line indicates the exact time 
delay To = 0.25, while horizontal dashed line indicates the level of filtering 
of singular values given by £ = 10 —13 . 


M 

To estimate 

M 

To estimate 

200 

0.40158 

700 

0.39113 

300 

0.25578 

800 

0.28392 

400 

0.19863 

900 

0.26543 

500 

0.14311 

1000 

0.20293 

600 

0.45358 

1500 

0.32837 


TABLE III 

Approximations of To in the four-pole example with 

PERTURBATION 10 -5 sin(107nr) USING EXTRAPOLATIONS WITH 
ORIGINAL FITTING REGIONS FOR VARIOUS M. THE EXACT VALUE 
To = 0.25, AVERAGED VALUE = 0.26586. 


made more focused around Tq by narrowing down the fitted 
region. Results of this procedure are shown in Fig. El This 
improves a little an average time delay to r p^ aver '> = 0.24216. 



Fig. 9. HEflHoo in the four-pole example with the added perturbation 
10“ 5 sin(107rx) and M = 200, 400, 600 and 800 together with their 
extrapolated quadratic fits constructed using more narrow fitting region. 
Vertical dashed line indicates the exact time delay To = 0.25, while 
horizontal dashed line indicates the level of filtering of singular values given 
In £ 10 ;:i . 

Next we show results when a smaller noise of amplitude 
a = 10~ 8 is added. The evolution of ||.Er||oo as T increases 
is shown in Fig. [10] We can see that the plateau error region 
in this case is at about 10 9 level, so the error growth 
region is bigger than in the previous case, which should make 
fitting and extrapolation more accurate. Indeed, extrapolated 
quadratic curves intersect the horizontal line with value f 
in a more localized region about Tq as shown in Fig. [TTl 



Fig. 10. Evolution of ||.Er||oo in the four-pole example with the added per¬ 
turbation 10 -8 sin( I Ott:/:) of smaller amplitude. The dashed line corresponds 
to the exact delay To = 0.25. 

while averaging of obtained approximation to To produces 
rp(aver) _ q 25436, which is more accurate than in the case 
with a higher amplitude a = 10 -5 . 



Fig. 11. ||T?r||oo in the four-pole example with the added perturbation 

10 -8 sin(107ra;) and M = 200, 400, 600 and 800 together with their 
extrapolated quadratic fits. Vertical dashed line indicates the exact time delay 
To = 0.25. while horizontal dashed line indicates the level of filtering of 
singular values given by f; = 10 1 


M 

To estimate 

M 

To estimate 

200 

0.46392 

700 

0.25502 

300 

0.27635 

800 

0.26081 

400 

0.25683 

900 

0.2444 

500 

0.26798 

1000 

0.25582 

600 

0.26391 

1500 

0.25292 


TABLE IV 

Approximations of To in the four-pole example with 

PERTURBATION 10 -8 sin(107ra;) USING EXTRAPOLATIONS FOR VARIOUS 

M. The exact value To = 0.25, averaged value T^ aver) = 0.25436. 


B. Transmission Line Example 

We consider a uniform transmission line segment with the 
following per-unit-length parameters: L = 7.574 nFl/inch, 
C = 2.61166 pF/inch, R = 16.278 mfl/inch, G = 5.58 
//S/inch and length C = 5 inches. The frequency is sampled 
on the interval (0,5.0] GHz. The scattering matrix of the 
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structure is computed using Matlab function rlgc2s. We 
consider the element H(w) = Sn(w) and impose the time 
delay To = 1.25 ns by multiplying H(w) by exp(— iwTo) to 
get the delayed transfer function H(w) = exp(— iwTo)H(w). 
The real and imaginary parts of H(w) are given in Fig. [TJ] The 



Fig. 12. HeH(w) and lmH(w) in the transmission line example. 

error curves for different M are shown in Fig. [13] indicating 
that the reconstruction error decreases quickly with M and 
reaches the level close to machine precision at M = 600. 
Constructing fitted quadratic error curves and finding their 




Fig. 14. Estimation of the delay time in the transmission line example using 
critical transition times T c as M varies. The dashed line corresponds to the 
exact delay To = 1.25 ns. 



Fig. 15. H-ErHoo in the transmission line example with M = 200, 400, 
600 and 800 together with their quadratic fits. Vertical dashed line indicates 
the exact time delay To = 1.25 ns, while horizontal dashed line indicates the 
level of filtering of singular values given by £ = 10 —13 . 

filtering affects the results more, extrapolation becomes less 
accurate. Averaging obtained approximations of T 0 produces 
rp{aver) _ ^ .2498 ns, that is very close to the exact value 
T 0 = 1.25 ns. 


M 

To estimate (in ns) 

M 

To estimate (in ns) 

200 

1.5194 

700 

1.2531 

300 

1.3147 

800 

1.2512 

400 

1.2793 

900 

1.2678 

500 

1.2506 

1000 

1.2348 

600 

1.2668 

1500 

1.2242 


Fig. 13. Evolution of ||2?fd|oo in the transmission line example as M varies. 
Vertical dashed line indicates the time delay Tq = 1.25 ns. 


intersections with jf-T^Hoo at T = 0 or finding times when 
these fitted error curves reach the value £ of the error for 
M > 600, we get a sequence of critical transition times T c , 
that we show in Fig. M Clearly, critical times T c converge 
to To and provide a good approximation of To for M > 500. 
Using an alternative approach when we extrapolate the fitted 
quadratic error curves to find their intersections with the error 
threshold £, we find approximations of To. Some of these 
curves for M = 200, 400, 600 and 800 are depicted in Fig. 
[T5l Approximations of To using extrapolation procedure for 
various values of M ranging from M = 200 to 1500 are given 
in Table [V] A good approximation of T 0 is obtained even with 
M = 300. As before, approximations of T 0 become better as 
M increases, but for very large values of M > 1000 when the 
reconstruction error falls below the filtering threshold £ and 


TABLE V 

Approximations of T 0 in the transmission line example using 
EXTRAPOLATIONS FOR VARIOUS M. THE EXACT VALUE Tq = 1.25 NS, 


AVERAGED VALUE TX 


= 1.2498 ns. 


C. Dawson’s Integral Example 

We consider here another analytic example |23l modeled 
by the transfer function 

H(w) = e~ iwT ° H(w ) 

where 

H(w) = e - ™ 2 — ^LD(w), 

V* 

D(w) is Dawson’s integral 

D(w) = e~ w f e* dt = e~ w erfi(u;) 

Jn 2 
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w 


Fig. 16. R eH(w) and lmH(w) in the Dawson’s integral example using 
N = 300 sample points. 



Fig. 17. ||.E fl ||oo in the Dawson’s integral example as M varies. Vertical 
dashed line indicates the time delay To = 0.125. 


and erfi(u') is the imaginary error function. Since 


H[ReH] =H(e~ w2 ) 



— ImH, 


function H[w ) is causal. Hence, the function H(w) is a causal 
function delayed with offset To. We use To = 0.125 and 
sample H(w) on [0,w max ] with w ma x = 20 using various 
numbers of points ranging from N = 100 to 600. Real and 
imaginary parts of H are shown in Fig. [16] The evolution 
of Halloo for various M is shown in Fig. [T7] where it is 
clear that critical transition times approach To. Constructing 
fitted quadratic error curves and extrapolating them to find 
their intersections with the horizontal line corresponding to 
the error value £ produces a set of approximation of To, 
shown in Table [VI] Averaging obtained approximations of To 
for M > 200, once some convergence is established, gives 

T (aver) = Q 12528 . 

It is interesting to note behavior of the relative error E r ^ 1 in 
this example. The evolution of its oo norm is shown in Fig. [18] 
It is clear that all profiles even for small values of M have 
a unique local maximum at T = Tq. 2-norm has a similar 


3 Please note that we use an opposite sign in the definition of the Hilbert 
transform than in 03 and EH. 


M 

To estimate 

M 

To estimate 

150 

0.16735 

240 

0.12362 

170 

0.13116 

300 

0.12661 

180 

0.12964 

400 

0.12233 

200 

0.12753 

500 

0.12578 

230 

0.12333 

600 

0.12414 


TABLE VI 

Approximations of T 0 in the Dawson’s integral example using 

EXTRAPOLATIONS OF FITTED ERROR CURVES FOR VARIOUS M. THE 
EXACT VALUE To = 0.125, AVERAGED VALUE FOR M RANGING FROM 
200 TO 600 is T^ aver) = 0.12528. 



Fig. 18. Evolution of the relative error H-E^Hoo in the Dawson’s integral 
example as M varies. Vertical dashed line indicates the time delay To = 
0.125. 


behavior. Even though the behavior of the relative error E^ 1 
can be used to determine the time delay in this example, we 
did not find the same pronounced behavior in other examples 
we considered. At the same time, extrapolating fitted quadratic 
curves of oo norms of the absolute error Er was a robust 
approach in all examples we considered. 


D. Stripline Example 

We simulated an asymmetric stripline modeled in lf37l with 
length L = 8 in, width IF = 14 mils, distances from the trace 
to reference planes Hi = 10 mils, H 2 = 20 mils, substrate 
dielectric Megtron6-1035, Laminate with a dielectric constant 
e r = 3.45 using a Cadence software tool with FEM full-wave 
field solver. The stripline was simulated on [0, w max ] with 
Wmax = 2 GHz. We analyzed element H(w) = Sn(w ) of the 
transfer matrix. The real and imaginary parts of H are shown 
in Fig. [T9] The evolution of ||.Er||oo for various M is depicted 
in Fig. [20] Even for high values of M, the error in causality 
does not go to machine precision and instead levels off around 
ICC 6 . This indicates that our finite element simulation results 
are accurate only within ICC 7 — ICC 6 . For causality char¬ 
acterization, this implies that data have noise/approximation 
errors with amplitude around ICC 7 — ICC 6 . Graphs of IjT^Hoo 
suggest that for M < 2000, the error is dominated by Fourier 
series approximation error, while for higher of M, the error 
is dominated by the noise/approximation errors from the finite 
element method. 

In this example, the time delay was estimated using a 
closed form microwave theory approximation as Tq = 8 x 
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Fig. 19. R eH and ImH in the stripline example with N = 1000. 



Fig. 20. Halloo in the stripline example for various M. Vertical dashed 
line indicates the closed form microwave theory time delay approximation 
To = 1.25809 ns. 


0.0254/(co/ v / e7) = 1.25809ns, where Cq = 3 x 10 8 m/s is the 
speed of light, 0.0254 is a conversion factor to convert from 
inches to meters. The error curves were fitted to quadratic 
curves as explained above. Because of relatively high errors 
in data, the fitted regions are not long enough. Besides, there is 
more nonlinear behavior of the error curves for high values of 
T > Tq . All this makes it difficult to estimate the time delay 
as shown in Fig. [2T] As can be seen, extrapolated quadratic 
curves do not focus at Tq but instead spread out around 
To similar to the four-pole example with an imposed noise 
considered in Section IV-AI This problem can be corrected 



Fig. 21. Extrapolated quadratic curves based on initial fitting range in the 
stripline example. 



Fig. 22. Extrapolated quadratic curves based on more narrow fitting range 
in the stripline example. The average time delay is T^ av ' r] = 1.2669 ns. 


by decreasing the fitting range and going more away from 
transition regions. The results are shown in Fig. |22] The 
approximations of To are given in Table IVIII Averaging them 
for values of M up to 3000 produces r j'( aver '> = 1.2669 ns, that 
agrees well with an analytically estimated time delay using a 
closed form formula. As in other examples, results with very 
high values of M > 3000, that are more affected by noise and 
approximation errors in data, are less accurate. 


M 

To estimate (in ns) 

M 

To estimate (in ns) 

80 

1.266 

700 

1.1987 

100 

1.2312 

800 

1.2179 

200 

1.2553 

900 

1.2593 

300 

1.2797 

1000 

1.2205 

400 

1.2826 

2000 

1.246 

500 

1.2413 

3000 

1.5878 

600 

1.1833 

4000 

1.0168 


TABLE VII 

Approximations of T 0 in the stripline example using 

EXTRAPOLATIONS OF FITTED QUADRATIC ERROR CURVES FOR VARIOUS 

M. The closed form approximation of the time delay is 
To = 1.25809 ns, averaged value for M ranging to up to 3000 is 

rp(av£r) = l2669ns _ 


VI. Conclusions 

We proposed a new method for time delay extraction from 
tabulated frequency responses. The approach uses the spec¬ 
trally accurate causality enforcement technique constructed us¬ 
ing SVD-based causal Fourier continuations, that was recently 
developed by the authors. The time delay is incorporated 
to the causality characterization approach by introducing a 
linear varying phase factor to the system of equations that 
defines Fourier coefficients. Varying time until a threshold 
time, that depends on the maximum frequency at which the 
transfer function is available, results in the reconstruction error 
between the given data and their causal Fourier continuations 
to go from an almost constant small value to a rapidly growing 
function at some critical transition time. The critical transition 
times depend on resolution and approach the time delay as 
resolution increases. Several sets of frequency responses with 
increasing resolution can be used to establish convergence and 
get an approximation of the time delay. Alternatively, when 
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only a limited number of samples is available, a growing 
portion of the error curve can be extrapolated to find an 
approximation of the time delay. The method is applicable to 
data that have noise or other approximation errors. A few sets 
of frequency responses can be used to improve the accuracy 
of time delay approximation by averaging the obtained results. 
The technique can be extended for multi-port and mixed mode 
networks. The performance of the method is demonstrated 
using several analytic and simulated examples, including data 
with noise, for which time delay is known exactly or can be 
estimated using other approaches. 
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Appendix 


column of V; K denotes the number of singular values that 
are discarded, i.e. the number of j for which rr t < £, where 
£ is the cut-off tolerance. 

It can be seen that constants Ai, A 2 and K depend only on 
the continuation parameters N, AI, b and £ as well as location 
of discrete points x 3 , and not on the function H. 

For brevity, we can write the above error estimate as 

H-ff — C(H + £)||_L 2 (f2) < €F + En + £T• 

Here 

e F = (1 + A 2y /N(M-K))\\H - H M lU^n) 

is the error due to a causal Fourier series approximation and 
it decays at least as fast as 0(M~ k+1 ), k is the smoothness 
order of H(x), which can be estimated numerically using 
reconstruction errors with different resolution (see 0). 

e t = Ai\/ AY6||iTM||L 00 (n=)) 

that is the error due to the truncation of singular values. It is 
typically small and close to the cut-off value £. 


Error Analysis of Causality Characterization 
Method Based on Causal Fourier Continuations 


In this section, we provide an upper bound of the recon¬ 
struction error between a given transfer function H(x) and its 
causal Fourier continuation C(H){x) in the presence of noise 
e in data. 

Denote by H m any function of the form 

M 

H m {x) = ^2,a k (t> k (x) ( 19 ) 

k =1 

where <j> k {x) = k = 1,..., M. 

Let A = UT.V* be the full SVD decomposition f34l of 
the matrix A with entries A k j = <p k (xj), j = 1,..., N, k = 
1,..., M, where U is an N x N unitary matrix, £ is an N x 
AI diagonal matrix of singular values a 3 , j = 1,... ,p, p = 
min(AT, AI), V is an AI x AI unitary matrix with entries V k j, 
and V* denotes the complex conjugate transpose of V. 

The following result is true 0. 

Theorem : Consider a rescaled transfer function H(x) de¬ 
fined by symmetry on 11 = [—0.5, —a] U [a, 0.5], where 
a = 0.5 , whose values are available at points x 3 £ 12, 

j = 1 ,... ,N. Then the error in approximation of II (x), that 
is known with some error e, by its causal Fourier continuation 
C(H)(x) defined in (0 on a wider domain fl c = [—b/ 2, b/2], 
b > 1, has the upper bound 


\\H - C(H + c)11 l 2 (0) < (1 + A 2y /N(M-K)) 

x \H — HmWl^iq.) + Hellion) )+A 1 ^Kjb\\H M \\L^) 

and holds for all functions of the form (IT9t . Here 


Ai 


.max | \vj(x) ||z, 2( n), A 2 


max 

3- Cj>£. 


11' v j( x )\\L 2 (n) 


and functions Vj(x) = Vkj(t>k{%) are each an up to M 

term causal Fourier series with coefficients given by the jth 


— (1 + A 2 \/N{M - Ar))||e|| ioo (n) 

is the error due to noise e in data. Numerical experiments 
reveal that e n has the order of noise e in data. 
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